




“The Jupyter Notebook is an open-source web application that allows you to create and share documents that contain live code, equations, visualizations and narrative text.” – jupyter.org

import sys
v = sys.version_info
print(f"We are using Python {v.major}.{v.minor}.{v.micro}")
We are using Python 3.9.5
GIS714: geospatial computation and simulation: graduate level course required for Geospatial Analytics PhD students

You can also try it out in Binder! (also linked from the GRASS GIS GitHub README)
# Import Python standard library and IPython packages we need.
import subprocess
import sys
# Ask GRASS GIS where its Python packages are.
sys.path.append(
subprocess.check_output(["grass80", "--config", "python_path"], text=True, shell=True).strip()
)
# Import the GRASS GIS packages we need.
import grass.script as gs
import grass.jupyter as gj
# Start GRASS Session
session = gj.init("data", "nc_spm_08_grass7", "PERMANENT")
# Make a new mapset for this assignment
gs.run_command("g.mapset", mapset="HW3_water_simulation", location="nc_spm_08_grass7", flags="c")
# Set region to the high resolution study area
gs.run_command("g.region", region="rural_1m")
And then we are ready to begin using GRASS!
Since grass.jupyter is written in Python, we'll use the GRASS Python API.
GIS714: Geospatial Computation and Simulation, spring 2022,
We'll use an example from the course on water simulation...
We create a map of flooded area with r.lake (GRASS manual for r.lake) by providing a water level and a seed point:
gs.run_command("r.lake", elevation="elev_lid792_1m", water_level=113.5, lake="flood1", coordinates="638728,220278")
# See results
flood1 = gj.Map(use_region=True)
flood1.d_rast(map="elev_lid792_1m")
flood1.d_rast(map="flood1")
# Display map
flood1.show()
Increase water level to 113.7m and 114.0m and create flooded area maps at these two levels.
#### Your Answer Here
# 3D map example
Calling display modules using Python magic short-cut:
d.rast -> Map.d_rast(), Map3D.d_rast
We will use two GRASS addons, r.stream.distance and r.lake.series, to estimate inundation with Height Above Nearest Drainage methodology (A.D. Nobre, 2011). We need to install them first:
#!g.extension r.stream.distance
#!g.extension r.lake.series
For this section, we will change our computation region to elevation which is a larger study area than we used above. Because we are using a new region (and also a higher threshold of 100,000), we need to run r.watershed again to compute the flow accumulation, drainage and streams. We convert the streams to vector for better visualization.
gs.run_command("g.region", raster="elevation")
gs.run_command("r.watershed", elevation="elevation", accumulation="flowacc", drainage="drainage", stream="streams_100k", threshold=100000)
gs.run_command("r.to.vect", input="streams_100k", output="streams_100k", type="line")
Now we use r.stream.distance with output parameter difference to compute new raster where each cell is the elevation difference between the cell and the the cell on the stream where the cell drains. This is our HAND terrain model.
gs.run_command("r.stream.distance", stream_rast="streams_100k", direction="drainage", elevation="elevation", method="downstream", difference="above_stream")
With r.lake.series, we can create a series of inundation maps with rising water levels. r.lake.series creates a space-time dataset. We can use temporal modules to further work with the data...
gs.run_command("r.lake.series", elevation="above_stream", start_water_level=0, end_water_level=5, water_level_step=0.5,
output="inundation", seed_raster="streams")
... or visualize the flood with TimeSeriesMap.
flood_map = gj.TimeSeriesMap(use_region=True)
flood_map.add_raster_series("inundation")
flood_map.d_legend(color="black", at=(10,40,2,6)) #Add legend
flood_map.d_barscale()
flood.show()
folium...

Moving data from GRASS GIS location projection to folium is a challenge
folium is projected in Web Mercator (EPSG:3857)
However, any coordinates (i.e. vector data) need to be specified in degrees of latitude and longitude (WGS84, EPSG:4326)
folium reprojects latitude and longitude to Web Mercator internally
We pass data to folium by reprojecting to temporary locations

import folium
# Create figure
fig = folium.Figure(width=600, height=400)
# Create a map to add to the figure later
m = folium.Map(tiles="Stamen Terrain", location=[35.761168,-78.668271], zoom_start=13)
# Create and add elevation layer to map
gj.Raster("elevation", opacity=0.5).add_to(m)
# Do some cool folium stuff!
# Like make a tooltip
tooltip = "Click me!"
# and add a marker
folium.Marker(
[35.781608,-78.675800], popup="<i>Center For Geospatial Analytics</i>", tooltip=tooltip
).add_to(m)
# and a circle
folium.Circle(
radius=120,
location=[35.769781,-78.663160],
popup="Great Picnic Area",
color="crimson",
fill=False,
tooltip=tooltip
).add_to(m)
# Add the map to the figure
fig.add_child(m)
# Display figure
fig
The InteractiveMap class allows users unfamiliar with folium to produce maps easily.
# Create Interactive Map
fig = gj.InteractiveMap(width = 600)
# Add raster, vector and layer control to map
fig.add_raster("elevation", title="Elevation Raster", opacity=0.8)
#fig.add_vector("roadsmajor")
fig.add_layer_control(position = "bottomright")
# Display map
fig.show()
Notebook format was generally well received:
There were several areas of the notebook deployment that needed improvement or that posed a difficulty to the class:
FOSS4G '22 GRASS GIS Workshop Materials (Anna Petrasova)
GIS714: Advanced Geospatial Computation and Simulation Course Repository
Helena Mitasova
Stefan Blumentrath
Vero Andreo
GIS714 Class Spring 2021
… and the GRASS community, NC State Center for Geospatial Analytics, Google Summer of Code